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Abstract 

We implement an efficient method of computation of two dimensional Fourier-type 
integrals based on approximation of the integrand by Gaussian radial basis func¬ 
tions, which constitute a standard tool in approximation theory. As a result, we 
obtain a rapidly converging series expansion for the integrals, allowing for their 
accurate calculation. We apply this idea to the evaluation of diffraction integrals, 
used for the computation of the through-focus characteristics of an optical system. 
We implement this method and compare its performance in terms of complex¬ 
ity, accuracy and execution time with several alternative approaches, especially 
with the extended Nijboer-Zernike theory, which is also outlined in the text for the 
reader’s convenience. The proposed method yields a reliable and fast scheme for 
simultaneous evaluation of such kind of integrals for several values of the defocus 
parameter, as required in the characterization of the through-focus optics. 
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1. Introduction 


The importance of the 2D Fourier transform in mathematical imaging and vi¬ 
sion is difficult to overestimate. For a function g given on in polar coordinates 
(p, 6) it reduces to calculating the integrals of the form 



( 1 ) 


in which the Fourier transform is expressed in polar coordinates {r,(p). For in¬ 
stance, the impulse response of an optical system (referred to as the point-spread 
function (PSF)) can be defined in ferms of diffraction infegrals of fhe form Q. The 
PSF uniquely defines a linear opfical sysfem and if is usually calculafed for a single 
value of fhe focus parameter ll2^ . 

Calculafing Ihrough-focus characferisfics of an opfical sysfem has a wide va- 
riefy of imporfanf applicafions including phase-diverse phase refrieval |[T3l IT4i . 
wavefronl sensing i44l . aberrafion refrieval in lifhography, microscopy, and ex- 
freme ulfraviolef lighf optics |[T5l[T^ . as well as in physiological opfics, where such 
calculafion can be used fo assess the efficacy of infraocular lenses 0131] im, sfudy 
fhe deplh-of-focus of fhe human eye 0|48l or defermine opfimal pupil size in reti¬ 
nal imaging insfrumenfs such as fhe confocal scanning laser ophfhalmoscope iflTl . 

Anofher rapidly developing field in imaging is fhe digifal holography ll25l |30l 
[40l@2l^ which finds applicafions in fhe quanfifafive visualization of phase objecfs 
such as living cells using microscopic objecfives in a digifal holographic sef-up. 
This technology of acquiring and processing holographic measuremenf dafa usu¬ 
ally is comprised of fwo steps, fhe recording of an interference paffern produced 
by a real objecf on a CCD, followed by fhe numerical reconsfrucfion of fhe holo¬ 
gram by simulafing fhe back-propagation of fhe image. Fourier-fype infegrals are 
an essential pari of several algorifhms used in digifal holography, for insfance when 
calculafing the Fresnel Transforms. 
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In this context, the numerical algorithms used to calculate the 2D Fourier trans¬ 
form play a central role. Unfortunately, such integrals can be explicitly evaluated 
in terms of standard special functions only in a limited number of situations. On 
the other hand, purely numerical procedures for these evaluations are computa¬ 
tionally cumbersome and require substantial computational resources. Thus, the 
most popular alternative are the semi-analytic methods that use an approximation 
of the integrand with functions from a certain class, for which some forms of closed 
expressions are available. This is the paradigm of the extended Nijboer-Zernike 
(ENZ) approach, explained briefly in Sectionj^ The ENZ theory constitutes an im¬ 
portant step forward in comparison to the direct quadrature integration. Its initial 
implementation ||8l 123 [27l was for circular pupils and for symmetric aberrations 
containing only the even terms of Zernike polynomial expansion (i.e., only cosine 
dependence), and presented some limitations such as its poor performance for large 
numerical apertures. Many of these shortcomings were fixed with the introduction 
of the vectorial ENZ theory for high numerical apertures ||9l, and in more recent 
publications ll45ll^l^ . 

In this paper we develop an alternative technique for calculating integrals of 
the form ([T]l by combining the essence of the ENZ approach with the approxi¬ 
mation power of the radial basis functions (RBE), which are standard tools in the 
approximation theory and numerical analysis, appearing in applications ranging 
from engineering and artificial intelligence to solutions of partial differential equa¬ 
tions, see for example iflSll . The Gaussian radial basis functions (GRBE) have been 
also used in the context of ophthalmic optics ll32l 13311341 . 

The structure of the paper is as follows. In Sectionj^we outline some facts from 
the approximation theory by RBE. The exposition here is very sketchy, and its main 
purpose is to mention some literature for the reader’s convenience. Sectionj^is de¬ 
voted to the derivation of the main formulas that can be used for the calculation of 
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general Fourier integrals ([T]l. The detailed algorithm is developed in Sectionj^for 
the diffraction integrals defined therein, including its efficient implementation and 
error estimates. An assessment of accuracy and efficiency of fhis procedure is car¬ 
ried ouf in Secfion]^ where if is compared also wifh fhe extended Nijboer-Zemike 
fheory (ouflined in Secfion]^ and wifh fhe sfandard 2D fasl Fourier Iransform. 

Some of fhe ideas underlying fhis approach have been presented in |[37]| . In 
this work we develop further the algorithmic aspects of this procedure, prove its 
convergence, discuss its efficient implementation, and compare its performance 
with some standard alternatives. 

2. Approximation by Radial Basis Functions 

Every semi-analytic method for calculating integrals is comprised of two phases: 
first we must approximate the integrand (or some of its factors) by a function from 
a given class, and after that the integral of the approximant is calculated in a closed 
form. This is the case of the ENZ approach (where the pupil function is approxi¬ 
mated by linear combinations of Zemike polynomials, see the details in Sectionj^, 
as well of the method put forward in this paper, which is based on approximation 
by a linear combination of Gaussian functions. 

A Gaussian radial basis function (or GRBE), 

is a strictly positive definite radial basis function (RBE) in with the shape 
parameter A > 0 and center at the origin. As a consequence, any set of data 
{ak,bk,Zk), with Pk = {ak,bk) G pairwise distinct, is unisolvent for Eagrange 
interpolation by a linear combination of shifts ^{x — ak,y — bk), that we briefly call 
Gaussian inferpolanf. 
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The theory of approximation by RBF (both the theoretical and computational 
aspects) has been rather well developed in the last 30 years, motivated by its ap¬ 
plications in numerical analysis (mesh-free methods) and neural networks, see 

e.g. MM- 

The accuracy of the Gaussian interpolants and the stability of the corresponding 
linear system depend on the number of data points and on the shape parameter A. 
For a fixed A, as the number of data points increases, the Gaussian interpolant 
converges to the underlying (sufficiently smooth) function being interpolated at 
a “super-spectral” rate ^), where h is the measure of a “typical” distance 

between the data points, and c is a positive constant, depending on A, see ||^|4^ . 
The selection of the appropriate A is important: smaller values of A yield better 
accuracy of the interpolant but affect significantly the stability of the linear system 
(the “trade-off principle”). Nonetheless, methods have been described to handle 
this situation lfT^l28l . 

Least-squares approximation by RBF has been developed in parallel within the 
artificial infelligence communify. For insfance, fhe supporf vecfor (SV) machine 
is a fype of learning machine, based on stafisfical learning theory, which contains 
RBF networks as special cases. In the RBF case, the SV algorithm automatically 
determines centers, weights, and threshold that minimize an upper bound on the 
expected test error BTl l43l . This has connections with other successful approxi¬ 
mation schemes, such as the adaptive least squares (see e.g. ifTSl Chapter 21], and 
also l[33l), the moving least squares ifT^ Chapters 21-27], and others. 

The discrete least squares approximation to a smooth function is known to 
converge with at least a power rate; for some error bounds and implementation 
issues we refer the reader to lITOl Chapter 8] (see also ifTSl Chapter 20]). 
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3. Fourier integrals for Gaussian RBF 


Consider the case when g in 0 is a Gaussian RBF with the shape parame¬ 
ter A > 0 and center at a point with cartesian coordinates {a,b), multiplied by a 
radially-symmetric complex-valued function /j; in other words, let 

x = pcos{d), y = psin(0). 

Using the polar coordinates for the center, 

a = qcos{a), b = qsm{a), 

we can rewrite it as 

g(p, 0) = (2) 

In this section we discuss an expression for o with g as in 0: 

1 poo r27i . . 

F{r,^-q,a)=- / +P ( 3 ) 

71 Jo Jo 

Theorem 1. For F defined in 0, 

F{r,^-,q,a) = 

where 

Q. = Q.{r,(^\q,a) = X^q^ + 27lirXqcos{^ — a) — Tl'^r^, (5) 

Iq is the modified Bessel function of the first kind, and denotes the Laplace 
transform. 

Proof. Notice that by 0, 

1 poo o r27t 

F{r,(j)-q,a) = -e-^‘^^ dph{p)pe-^P / g2Ap^cos(e-a)^2;r;rpcos(e-0)^p_ 

71 .Jo .Jo 

( 6 ) 


loilV^hi^-) (A), 


(4) 
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For arbitrary constants A,B G C we have the following identity: 

^2Acos{e-a)^2Bcos(e)^Q ^ (^2y/A2 + 2ABcosa + B2^ . (7) 

Observe that Iq is an entire even function, so that the value in the right hand side of 
Q is independent of the choice of the branch of the square root. 

Identity (|7]l is a straightforward consequence of the fact that 

Acos (6 — a) +Bcos(6) = \/a2”-|-^AScos^H-^cos(0 — 7 ), 


where 


7 = arctan 


A sin a 


^B+Asina J ’ 

(with an appropriate choice of the branch of the square root), along with the identity 
(|T|, formula (9.6.16)): 


/o(z) = l 

TT 7o 


Using 0 in Q we conclude that 


F{r,<p-,q,a) = 2e y e Iq h{p)pdp 


l^e-^Plo(2y^y{^)dp, 


( 8 ) 


which is equivalent to 0 . 


Formulas in Theorem[T]are the basic building blocks for the algorithm proposed 
below. Let us particularize these identities to the case of a circular exit pupil of an 
optical system, when h can be taken as the characteristic function of the interval 
[ 0 , 1 ]: 


Corollary 2. Let 


Kp) = X[0,i]ip) = 


1 , if0<p<l, 
0 , otherwise. 


1 






(9) 


Then F in Q takes the form 


F{r,^-,q,a)=e ^(r,0;^,a)^ 

s=0 


where 


/■i 


ms{X) = < 


( 10 ) 


/ e-^Pp^dp, A/0, 

Jo 

(1+5) \ A = 0 . 

Furthermore, \ms{Ji) \ < 1/or Re(A) > 0, and the series in Q is uniformly conver¬ 
gent for all values ofr,q>0, (/>, a G M and A > 0. 

Remark 1. Integration by parts and straightforward calculations show that for s = 
0 ,1,2,... and A / 0, 


1-0 ^ {s+l)m,{X)-e ^ 

mo(A) = —r-, m^+i(A) =- 


A 


( 11 ) 


Notice also that 


( d \ 

~lx) '”o(A), s>l. (12) 
Proof. By Q, the crucial step is the evaluation of the integral 


e-^Ph{l^p^ dp. 


to 


The Taylor series for the Bessel function (see d, formula (9.6.12)), 

ieo(^ 

is locally uniformly convergent on the whole plane, and thus 


/o(2 


= i 7mP' 


/ \ °° 
e-^Plo (2 dp = Y, 


{m,{X)Tiy 


to t-y 


(13) 


Furthermore, the bound for m^X) is trivial for purely imaginary A, while for 
Re(A) >0, 

l_g-Re(A) 


\ms{X)\ < |m,(Re(A))| < f\-^<^^Pdp = |mo(Re(A))| = 

Jo 


Re(A) 


< 1 , 


which concludes the proof. 









Clearly, formulas of Theorem [T] and Corollary can be easily extended to the 
case when g is a linear combination of functions of the form Q, that is, 

g{p,e) = hip) f ^ 

k=\ 

where c^: G C are certain coefficients, Xk>0 are the shape parameters, and {q^, CCk) 
are the polar coordinates of the corresponding centers of the Gaussian RBFs. The 
corresponding Fourier transform F, defined in ([T]), can be wriffen as 


K 

I 

k=i 


Fir,(^) = Y,(^lce ^ ^irA'^qk,CCkY, 


to tf 

wifh Q. defined in Q. 

If we assume additionally fhaf all fhe shape paramefers are equal, Xi = ■■■ = 
Aif = A, we can rearrange fhe calculations as follows: 

fisir,^) = Y,Cke^^FQ,{r,(l)-,qk,ak)'. 

to t) ti 


4. Calculation of diffraction integrals 


Optical systems can be described by means of the complex-valued pupil func¬ 
tion Pip, 6) that, in the case of a circular pupil, can be expressed in (normalized) 
polar coordinates as 

Pip,e) =A(p,e)exp(^-/^iT(p,e)^, (i4) 

where A(p,6) is the aperture or amplitude transmittance function of the optical 
system, VF(p,6) is the wavefront error, and constants n and £ denoting the re¬ 
fractive index and the wavelength of the light, respectively. According to Fourier 
optics fTl, the complex-valued point-spread function of such a system is given by 
the diffraction integral 

Uir,^-,f) = - / F(p,/)P(p,e)exp(27r/prcos(e-(/))) pr/ejp, (15) 

K Jo Jo 
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where / represents the defocussing {f = 7i/2 corresponds to one focal depth), and 
(r, 0) denote the polar coordinates in the image plane. The focal factor F is 


^(P,/) =exp 



(16) 


V / 

where 0 < < 1 is the numerical aperture (NA) of the optical system. For small 

values of (say, < 0.5) the focal factor is well approximated by exp(//p^), and 
the resulting integral is 

U{r,^-,f) = — / exp{ifp^)P{p,d)exp(27iiprcos{d — ^))pdddp, (17) 

71 Jo Jo 


known as the Debye approximation of the Rayleigh integral of the optical impulse 
response or the point-spread function (PSF) of the system, evaluated by 


PSF{r,(t,-f) = \U{r,^-f)\\ 


In this section we discuss the computational framework and the efficient im¬ 
plementation of the algorithm for calculation of such diffraction integrals for small 
NA, using ( fTT] ), although our method easily carries over to any NA and to the focal 
factor ( [T6l ). 

Observe that formally the defocussing term (exp(//p^)) in ( fT7| ) could be ab¬ 
sorbed in the aberration term (pupil function). However, from the physical point of 
view, the aberration term is an intrinsic error of the optical system, while the de¬ 
focussing is a deliberately introduced defect, which can take on several values that 
are relatively large with respect to the aberrations of the system. For this reason, a 
separate treatment is commonly preferred, see O. 


4.1. Computational framework 

In practice, the wavefront is sampled at a discrete and finite set of points on the 
pupil (the unit disk). In other words, the input data set has the form {xj,yj,Wj), 
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with Wj = W{xj,yj). Values Wj can be measured directly or obtained by standard 
procedures. 

As it was mentioned above, semi-analytic methods for calculation of Fourier 
([T]) or diffraction integrals ( fTT] ) are based on an approximation of the integrand by a 
function of a suitable chosen form. The brief outline of the results in Section|2]and 
the reasonable assumption that the complex pupil function P in ( fid) ) is sufficiently 
smooth allows us replace the actual pupil function by its approximation by a linear 
combination of Gaussian radial basis functions (GRBF). It is convenient to point 
out that the implementation of this approximation and the achieved accuracy are 
not subject of this paper, and are treated in the literature, see Sectionj^ 

Hence, our starting point for calculation of ( fTT] ) is the ansatz that the complex 
pupil function P is given by an expression of the form 

p{p,e) = f c,g,(p,e), g,(p,e) = 

k=i 

where c*. G C are certain coefficients, > 0 are the shape parameters, and {q^, a^) 
are the polar coordinates of the centers of the Gaussian RBFs. A practical way 
to obtain ( [T^ is to fix a basis {gk}k=i of GRBF, or equivalently, to choose, for 
each index k, values for < 7 ^, and The optimal choice of the shape parameter 
based on the given data is a highly non-linear problem, usually solved by cross 
validation, see ll20l 1^ . Using this method we arrived at the following rule of 
thumb been tested in practice: take the same shape parameter for all k’s, 

Ai = ---=A^ = Ag [10,20], (19) 

select r G N and create a grid of r x r equally spaced points in the square [—1.2,1.2]^. 
The goal of covering an area larger than the unit disk is to deal with the Gibbs phe¬ 
nomenon at the boundary of the disk. The polar coordinates of these K = r^ points 
will be the pairs {q^, CCt) in ( [T^ . 
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With the basis {gj,} chosen, we can calculate the coefficients Ck fitting the pupil 
function by means of the linear least squares. For a high number of centers K this 
problem is ill-conditioned, and the use of Tikhonov regularization |[24l is recom¬ 
mended, with an appropriate choice of regularization parameter, for instance, by 
the L-curve method or by Morozov discrepancy principle. In our experiments, im¬ 
plemented in Matlab, we used the algorithms described in fTTl . see also |[29l . For 
the error estimates, see the references in Sectionj^ 

An alternative approach to the 2D nonlinear approximation of functions by 
sums of exponentials is described in O. 


With all the parameters in the representation ( [T8| ) calculated we can use the 
formulas obtained in Sectionnamely 


K 


U{rA'J) = £ ms{Xk if) 

k=l s=0 


( 20 ) 


with 


D.{r,<p-,q,a) = -|- 27 r/rA( 7 cos((/) — a) — 

Furthermore, with the additional assumption ( [T^ we can rewrite this formula as 

^ D(r, a-t)'- 

.to (^!r 

( 21 ) 


4.2. Efficient implementation of the calculations 

The recurrence relation o can present some numerical instability, especially 
for large value of the imaginary part of the argument; these considerations suggest 
to replace the coefficients mffX), defined in ( fTO] ), by fheir normalized counferparfs, 

m^(A) 


inffX) = 


( 5!)2 ’ 


fhaf can be efficienfly generafed for A / 0 by fhe following double recurrence: 


m,+i(A)= 


fhs{X) 

5-1-1 


-t,(A) 


Ts+l (■^) — 


t.(A) 

(5 + 2)2 


5 = 0,1,2,... (22) 
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with the initial conditions 


To{X) = e-\ = 


(23) 


Since formulas in (| 20 l)-([ 2 T]) contain an infinite sum, we must choose a cut-ojf 


parameter, S G N, so that the expression in ( |2()| ) is replaced by 

K. S 

- if)^{r,^ ,qk,aky, (24) 

k=\ i=0 

while ( |2T] ) becomes 

s 

U{r,(^-,f) = J^fhyXk-if)Hyr,(j)), (25) 

s=0 

where we denote 


= '^Cke ^'‘‘^kQ_{r,^,qk,aky. 

k=i 

It is worth observing that in these formulas the defocus parameter / and the coor¬ 
dinates r and (p are independent: / does not appear in fl, while m^Xk — if) need to 
be calculated only once regardless the number of points at which we evaluate the 
functions. 

In practice, we want to find the values of U{r,p-,f) at a vector of J points 
on the plane, given by their polar coordinates (r, 0 ), and for a vector of defocus 
parameters /. In other words, we need to compute efficiently the matrix 


U = {U{rj,Py,U))^.^C^-f 

Let us discuss the algorithm under the assumption ( [T^ , so we will use formula ( | 2 T| ). 
Applying ([22l)-([2^ to the vector /, we obtain the matrix 

M = (myx - //.))r=i':;i e 

which, as it was mentioned, is computed once for all {rj, (pj)'s. 
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Another observation that will speed up computations is that if P and Q are ' 
vectors of Cartesian coordinates of points on with polar coordinates (r, 0) and 
{q,(x), respectively, then by Q, 

a) = X^q^ + 27 lir?iqcos{^ — a) — Tl^r^ = {XQ + KiP)^{XQ + TliP), 

(26) 

where the superscript T denotes the transpose of a matrix. In other words, Q.{r, ^^q^a) 
is the square of the euclidean distance (in M^) from XQ to —niP. It allows us to 
encode efficiently the evaluation of = {0.kj) e 

^k,j — jt Of/;), 

by first finding fhe Carfesian coordinafes of fhe confers, 

ak = qkCOsak, bk = qt^inak, k=l,...,K. 


Finally, we find fhe column vector 

d = {cke-^^^k=i . ^^1) 


We describe fhe resf of fhe procedure in fhe Algorilhm[T] 

If should be noticed finally fhaf we can slighfly increase fhe accuracy adding a 
consfanf ferm fo fhe approximation of fhe complex pupil function, fhaf is, insfead 


of ( [181 ) using fhe represenfafion 
K 


F(p,0) = co + £c,g,(p,0), gk{p,e)=e-^^^^^+P -2p«cos(e-«,))^ 


k=\ 


where co can be faken as fhe average of fhe values of P af fhe given poinfs. No¬ 
tice fhaf fhe confribufion fo [/ corresponding fo fhis ferm can be evaluafed using 
Algorifhm[^wifh A = 0. 


14 



input : Matrices M G e and d G 

output; Matrix U 

Initialization: 

Set t7 = 0 G 

Define R G as the matrix of I’s, and compute 

H = d^R. (28) 


// At this stage, H = ji G 

for 5 ^ 0 to S — 1 do 

Calculate 

U ^U + M{:,s+l)H, 


Update R by 


R i — R. =t= ri 


// We use the Matlab notation (.*) for the 
term-by-term multiplication; 

Compute H using ( |2^ ; 


Algorithm 1: Implementation of formula (21 1 . 
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4.3. Error estimates and convergence analysis 

The algorithm described above presents two main sources of errors (besides 
the standard truncation and machine arithmetic errors that we do not discuss here): 


the fitting error in ([18]); 


the truncation error, given by the choice of S, in replacing the infinite series 
in ( |20| ) or ( |2T| ) by finite sums. 


The former is not subject of this paper, and the reader is referred to the discussion 
in Section |2j 

Regarding the the truncation error, observe that by Corollary]^ for A > 0, 


ms{X-if)\ < 1 . 


On the other hand, by expression (j^, 

| 0 | = \D.{r,^-,q,a) \ = \J _|_ (27trXqco&{^ — a))^ = X^q^ + Tl^r^. 


In consequence, the truncation error in (20 1 is given by 


/ , / ■',2 idk^OCk) 


< I 


1 


i=s+l 


(.!)2 




or in other words, the remainder of the Taylor expansion for the modified Bessel 
funclion 


lo ( 2\lX^ql + n'^r'^ 


and fhus. 


L - /o2 ^{r,(^,qk,Ctk) 

s=s+i 


< 


[X^ql + Tl^r^) 


.2„2\5+1 


max 


(S + 1)! 0<t<Xlql+-n34 


r(S+l) 


(0 


For fhe practical implemenfafion, in our algorifhm wifh ( [T^ , ^ < 1.7, r < 2, 
A < 20 , so fhaf 

|0| < 1200. (29) 


16 

















This bound is very conservative, but shows that in the worst case scenario S = 100 
provides a truncation error of order 10^^. 

Although, as the numerical experiments explained in Sectionj^show, the method 
performs well even for reasonably large values of r and / in ( [TT] ), we can slightly 
improve accuracy by using the scheme ( |20l ), and replacing the infinite series by its 
Fade approximant of a suitable order with center at the origin (instead of a mere 
truncation, as we did above), see e.g. 0. 

5, Outline of the extended Nijboer-Zernike theory 

In Section|^we are going to produce the results of some numerical experiments 
regarding the behavior of the algorithm described above, and to compare them with 
the standard procedures used for the computation of diffraction integrals. 

One of the best-known semi-analytic methods of calculation of these integrals 
is the so-called Extended Nijboer-Zernike (ENZ) theory l|8]|2^, similar in spirit to 
the method presented in this paper, with the main difference that the pupil function 
P is expanded in terms of Zemike polynomials, instead of Gaussian radial basis 
functions. 

Recall that the (n,m) Zernike polynomial is defined as 



N'”Rn (p)cos(m0), m>0 
A“/?lr'(p)sin(|m|0), m<0 


where 0 < p < 1 and 0 <6 <271 are polar coordinates, and fhe double index {n,m) 
is such fhaf n>0,m< \n\, and n — niis an even number. The radial parf (p) is 


a rescaled Jacobi polynomial 
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given explicitly by 



P 


n—2s 


(30) 


With the normalization constant = \/{2 — 5o^m){n + 1) the Zernike polynomi¬ 
als become orthonormal on the unit disk D = {(x,y) G 'M?/x^+y^ <1} with respect 
to the area measure (see e.g. n Section 9.2] for further properties of Zemike poly¬ 
nomials). 

For the reader’s convenience, in this section we present an outline of the ENZ 
formulas. Assume that (with a suitable approximation) 


/’(P,0) = L«(P,0)- 


(31) 


If the coefficients in the expression ( [3T| ) are known, then the integral U in ( [TT] ) 
takes the form 


[/(r,(/);/) = £OT(r,(/.;/), 


(32) 


where 



According to the original ENZ theory, for the double-index {n,m) with m > 0 (a 
necessary condition for the applicability of these formulas) we have 


Un (r, </>;/)= (r,/) cos(m(/)) 


(33) 


with 



(34) 


and 


T> / \ ^ Jm+k+2j+\{27lr) 

Bk{r)=l^Ukj -—-- 

;=0 


(35) 
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with the coefficients 


“‘' = <-‘>'TnT7Tri jl, 


q + k + j\ 
k ) 


(36) 


where p = {n — m)/2 and q = {n + m)/2. Here Jy is the Bessel function of the first 
kind and order v, see e.g. |[T] Chapter 9], and formulas (|34l)-([^ are known as the 
power-Bessel series (ENZ-PB) expression for P™. 

Observe that these formulas contain a removable singularity at r = 0, so some 
extra care when organizing the calculations should be put. Moreover, they also 
present difficulties when the defocus parameter / is large. According to O, ap¬ 
proximately 3I/I terms are needed in general to accurately evaluate expression 
( [3^ , and the practical use of that formula is limited to a range |/| < Sn. 

In order to overcome these difficulties an alternative approach has been put 
forward in ll27ll . The main idea is the use of Bauer’s identity: 


exp(//p2) = exp £ {2k+\)i^jk ^m(p)> (37) 

where is the radial part of the Zernike polynomials, and 

A(z)=(^) Jk+{i/ 2 )iz), k = 0,l,... (38) 

Using ( |3T] ) and ( [37| ) in ( [TT] ) we reduce the problem to the computation of in¬ 
tegrals containing products of two radial Zemike polynomials, and . The 
crucial idea is to use the linearization formulas for these products of the form 


nnli nm2 

‘^m\ +2pim2+2p2 


/ 


mi+m 2 

mi+m 2 + 2 /’ 


where coefficients ci defined as 


TMi ^m2 mi+m2 
2^ j Pl,SlJ P2,S2^Si+S2-2t,l 
suS27 
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are expressed in terms of the quantities fp\ and g'"i given explicitly by 


fp,s = i-^) 


p-s^ 


2s+\ 


p + s+\ 

for 5 = 0,... ,p, and 

^ m + 2 / + 1 
~ m + u + l + l 


m + p — s—\\ fm + p + s 
m—l / \ s 


p + s 
s 


m 

u — l 


u -\-1 

I 


m + l + u 
m + l 


for M = + m. For m = 0, expressions (391 and (401 boil down to 


(39) 


(40) 


fp,s = ^p,s, = 


where 5 is the Kronecker’s delta. 

As a result, we obtain the so-called Bessel-Bessel series (ENZ-BB) expression 
for F”’ in ( |34| ): for n, m nonnegative integers with n — m>0 and even, with p = 
^(n — m) and q = ^(n + m), 


VnirJ) = exp f hf \ £ {2k + 
/ k=0 


f ^ 'V ( iVu' 


(41) 


where Iq = max{0, k — q,p — k} and 

Wk,l = ^ ^ fp]s^k,s,tgk+s-lt,li 

s=0 t=0 

with 

_ 2^1 +2s2 — 4t + 1 /\ 

“ 2^i+2^2-2t + l V J ’ 

for t = 0,..., min{5i, 52 } and Ak = {^^). 

In the special case of m = 0 we have that 


(42) 


(43) 


Wk,k+P-2j = h,p,j, 7 = 0,l,...,min{k,p}, 
while all other Wkj vanish. In this way, all Wkj > 0 and Wk.i = 1 ■ 
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The conclusion of 11271 is that this new scheme is valid and accurate even for 
very large values of /. By analyticity, the expression ( |4T] ) remains valid also for 
complex values of /, which is important in a number of practical applications, 


with the downside in the removable singularity of ( |38] ) at / = 0. In practice, it is 
convenient to use both ENZ schemes if the defocus ranges from 0 to large values 
off. 

Further enhancement of the Bessel-Bessel approach was suggested recently in 
ISIESSTI, where equation ( |4T] ) is rearranged for its evaluation in the following 
way: 

oo n-\-2k 


k—0 h——n+2k, 

h—m even, 
h>\m\ 


(44) 


Functions C and B of the defocus parameter / and of the radius r, respectively, are 
given by: 


Cikif) = {2k + l)i exp -if jk -f 


1 


1 


Buir) = 


Jh+ijlKr) 

2nr 


(45) 

(46) 


and the coefficients are defined as 






y«2 

fm2 




(47) 


where fhe ferm in parenfhesis above represenfs fhe Wigner 3j-symbol (see 
Chapter 34]), and is a real number. These expressions consfifufe fhe “enhanced” 


Bessel-Bessel series (ENZ-EBB) expression for V”’ in (34i. 


6. Comparative assessment of accuracy and efficiency 

In this section, the performance of the methods based on the Gaussian radial 
basis functions (GRBF) proposed in this paper is analyzed and compared to the 
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popular alternative approaches. 

Since the closed analytic expression for a Fourier transform type integrals like 
([T]) and ( [TT] ) is possible only for most elementary integrands, for their computa¬ 
tion we must rely either on numerical or on semi-analytical methods (or analytical 
approximations). The first group comprises quadrature-type numerical procedures 
in which the integration over a 2D domain is replaced by the calculation of a dis¬ 
crete sum based on evaluation of the integrand at a discrete set of points, with the 
particular challenge of integrating a highly oscillatory function. The most popular 
approach is the bi-dimensional discrete Fourier transform (we will refer to it as the 
FFT2-based approach), calculated via the Fast Fourier Transform algorithm. Its 
efficiency can be subsfanfially enhanced using fhe fracfional Fourier fransform [!5l| 
or a “bufferfly diagram” ideas |[T2l . see also |[ 22 l . 

The besf known represenfafive of fhe second group is fhe Exfended Nijboer- 
Zernike (ENZ) fheory fhaf was briefly exposed in Section 

In fhis section, fhese alfernafives are discussed and compared wifh our mefhod, 
paying special affenlion fo fheir compufafional complexify, precision, accuracy and 
speed. We use fhe “naive” notion of complexify, undersfanding by fhis fhe number 
of real fioafing poinf operafions (flops) needed fo run fhe algorifhm. Since fhe exacf 
number of flops is in general difficulf or nol feasible fo calculafe, fhe leading ferm 
for large values of fhe parameters is used. 

The assessmenfs are performed under assumptions allowing fhe application of 
all fhe mefhods explained above. In ofher words, we will evaluate fhe diffracfion 
integrals U = U{r,(j)-,f), defined in ( [T7| ), in a grid of poinfs (in polar coordinates) 
(r, 0) G and for a vector f G of values of fhe defocus parameters /. 

Eurfhermore, for fhe ENZ formulas we fake info considerafion only fhe terms con¬ 
taining the cosine, so that the original ENZ-PB formulas can be used. 
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6.1. Some remarks about FFT V5 semi-analytic methods 

In the FFT2-based scheme, the value of U {r,^-,f) is computed by means of the 
bi-dimensional fast Fourier transform. The crucial step is the substitution of the 
double integral in ( fTT] ) by a discrete sum. Notice that each new value of the defocus 
parameter / obliges to calculate the values of U completely, at a computationally 
high cost. Furthermore, the use of the FFT requires re-sampling the wavefront at 
a regular Cartesian grid covering the pupil; for convenience, the length of the grid 
should be an integer power of 2 in each direction. 

Another remarkable issue to be addressed when using the FFT2 scheme is the 
aliasing, a typical phenomenon that can appear due to discontinuities of the inte¬ 
grand. In order to prevent this, the pupil must be small in comparison with the 
sampled area (or in other words, we must extend the pupil to a larger region, set¬ 
ting the pupil function to zero in the complementary domain), resulting in a large 
area where U{r,^-,f) is negligible Il22l . Thus, a big portion of the computational 
load of this scheme is useless, and in general the spatial resolution needed with 
this method will be much higher than that required for the semi-analytic methods 
like discussed here. However, FFT2 is the standard method used in commercial 
ray tracing packages, such as Zemax or Code V. 

The advantage of the semi-analytic approaches, such as the ENZ theory or the 
method proposed here, is that they reduce the computation of U in ( [T7| ) to evalua¬ 
tion of more or less complex explicit expressions in terms of some elementary or 
special functions. One of the benefits of having these formulas is a better control 
of the image domain being computed, increasing the precision. Another impor¬ 
tant advantage is the huge boost in performance gained when a parallelization of 
calculations is done for multiple values of the defocus parameter /. 
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6.2. Computational complexity 

Let us analyze and compare the computational cost of evaluating the diffraction 
integral using the ENZ theory and by the GRBF scheme, for a grid of values r G 
4> G (so that with the notation of Section 4.2 J = N'^), and for M distinct 
defocuses, gathered in a vector f G M^. 


The ENZ—PB method 


We start with the basic ENZ method, corresponding to formulas ([32|)-([36|). Ob¬ 
serve that in ( [^ we must compute values of the Bessel functions of the first kind 
Jy, usually by their power series 




(48) 


(see e.g. [formula (10.2.2)1 |[35l ). In practice this series must be truncated at some 
term B — \. For the values of variables and parameters usually appearing in this 
context, we have determined experimentally that we can ensure the truncation error 
of 10^^ taking B = 15. The computational complexity of computing the Gamma 
function at a value of Td (10') is G (t^log(f)^). With B = 15 and with the usual 
Zernike fit (up to the 8th order) it is sufficient to take f < 3, which yields a max¬ 
imum of 11 operations for the Gamma function evaluation. Thus, for N different 
values of r we need roughly ^(3N-|-25) operation to evaluate a single term in 


the series ( |48] ), and in consequence, the computational cost of evaluating a single 
Bessel function is G {3BN). 


With the same assumptions, the evaluation of each coefficient in ( |36| ) is at most 
^ (150) flops. Consequently, the complexity of computation of each B^ defined in 


IS 


^((p + l)(30B + 186)) 


operations. 
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In order to evaluate expression ( |^ , we truncate the infinite series at some term 
5 — 1 , in which case the number of operations required to evaluate by ( [^ is 

^(5((p + l)(30B + 186)+2)+5 + 8 ) = ^(30B5(p +1)). 

Propagating this estimate to ( [3^ , computed also at N different values of 0, we 
need about (30BS(p + 1) +1N'^M + IN) flops for its evaluation. 


In order to estimate the overall cost for expression ( |32| ) we have to take into 
account the indices n and m therein: n ranges from 0 to a certain n* (the maximum 
radial order), and for each fixed n, m takes even/odd non-negative integer values 
< n. Thus, at each level n, only G (njl) Zernike polynomials are used (recall the 
restriction of m > 0 for ENZ formulas). Then, p = — = G (n/4) in average. 

In addition, the maximum radial order n* yields a total of (n* + 1) Zemike 
polynomials, from which only those with the cosine are used, so that the effective 
number of polynomials is approximately K = ^n*{n* + \). 

Combining all the previous assumptions, the complexity of computing 17 at a 
grid of 77 X 77 X M points for (r, 0, /) by the ENZ formulas can be estimated to be 

cost{UENZ-PB) = G (^77^/^5775+ 77 (77^77 + 77775 + 5775)) , 


where 5 and 5 are the truncation values for the series in ( |48] ) and ( |34l ), respectively, 
and K is the total number of Zemike polynomials included in the pupil representa¬ 
tion ( [3T] ). 


The ENZ—BB method 

We can perform a similar analysis for the improved ENZ scheme, correspond¬ 
ing to equations ([^-([42l). The only difference with respect to the previous discus¬ 
sion lies in the computation of E™, so we only need to re-estimate the computa¬ 
tional cost of E™ in equation ([4T]). 
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With the same assumptions used before, the evaluation cost of coefficients in 


equations (|39|), (|4f)|), (|43|) is of at most of 110 flops. Thus, the complexity of 


coefficients ; in formula ( [4^ is of ^ (330p^) flops. Recalling the estimate of 
G {3BN) flops for each Bessel function Jy, the cost of the finite sum (that with 
index 1) in expression ( |4T| ) is G {{k + p){3BN + 330pk + 3N)). 

The cost of evaluation of the Bessel spherical function 7^ is G {3BM + 3M) = 
G {3BM). Thus, if the improved formula for is used with the truncation of the 
infinite series at the term 5 — 1, its complexity is 


G {S {3BM + BN + {S/2 + p) {3BN + 165pS)) + 2M). 


Assuming as before that p = ^{n — m) = G{n/4) in average, it can be simplified 
to 

G (^nS^ (215+ lOn) + l.5BNS^ + ^nSBN+ NMS+ BMS^ . 

Since the rest of the calculations is the same as in the basic ENZ scheme, the 
complexity of evaluating 17 in a grid ofNxNxM points (r, <p,f) by the improved 
ENZ formulas (|3^-(|4^ is 

COSt{UENZ-BB) =G + {S^ + bnS) + 

+ K {N^M + NMS + BNS + S^ + S'^BN) + , 


where B and S are the truncation values for the series in (481 and ([34|), re 


spectively, and K is the total number of Zernike polynomials included in the pupil 
representation ([^. 


The ENZ-EBB method 

This scheme is almost the same as in the ENZ-BB method, except for the rear¬ 
rangement of the computation of the term +"7 Thus, only the complexity estimate 
of this term needs to be updated. 
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With the main basic assumptions made before, and since the Gamma function 
is involved in the Wigner 3j-symbols, the number of operations needed to compute 
each coefficient A'^^\n2%7^ of about 250. Recalling the estimate of ^ {3BN) flops 
for each Bessel function Jy, the cost of the innermost sum in expression ( |4^ (the 
finite sum of index h) is {3nBN + 250n + 3A^)). 

As in the previous case, the cost of evaluation of the Bessel spherical function 


jk is (3BM) and then the complexity of expression (^i, when truncating at the 
term 5 — 1, is approximately 


ff{S{3BM + n{3BN + 26B + SN))+MN). 


Combining the cost of R™ with that of the rest of the ENZ BB scheme studied 
above, the complexity of evaluating 17 in a grid of N x N xM points by 

the enhanced ENZ BB formulas (|44l)-(|47|) is 

cost{UENZ-EBB) (BNS+NS + BS) + 

+ K {N^M + NMS + BNS + BMS) + N^m'^ , 


where B and S are the truncation values for the series in i and ( |34| ), respectively, 
and K is the total number of Zernike polynomials included in the pupil representa¬ 
tion ( [3T] ). 


The GRBF method 

Now we switch to the computational complexity of the GRBE based formulas 
described in Section Again, we are interested in evaluating expression ( |24| ) at 
a grid r G M^, (jj £ and / G M^, with the additional assumption ( fT^ , so that 
formula (|25|) is used. 


According to the discussion in Section 4.2 once the cut-off parameter S has 


been chosen, the double recurrence (22 1 -( 231 must be evaluated for s < S, with an 
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estimated computational cost of {ASM) flops. After that, the calculation of Q. in 
( [2^ takes about {IN^K) flops. The computational load of finding d by ( [Z/] ) is 
almost negligible, requiring only {AK) flops. Recall that Q. and d are computed 
only once during the initialization of the algorithm. 

The computational complexity of evaluating each term Hs in the mesh is of 
^ {K + 2N^) operations, yielding a total of {KS + 2SN^) flops for the evaluation 
of all the required Hg terms. 

Summarizing, the cost of calculating U in (|24l) according to Algorithm 1 is of 


{2N^MS + TKN^ + N^M) . 


Recall also that the minimal estimated computational complexity of the FFT2 
method for a single value of /, even for optimal implementation, is of (A^ log(A)), 
which corresponds to the cost of the FFT, the most computationally demanding 
part. 

Table[T]summarizes the leading terms of the computational cost of each method, 
with the assumption that the same values of B and S for all semi-analytical meth¬ 
ods have been chosen. A quick comparison shows that the GRBF approach is much 
more efficient than the ENZ theory, especially for a large amount of values for /. 
The FFT2-based method seems to be of similar complexity to GRBF, but in prac¬ 
tice the number of sample points N needed to achieve a reasonable accuracy for 
FFT2 will be much larger than that required for GRBF, which in our experiments 
was set to 400 (see the description in the next section). 


6.3. Execution time 

Since the complexity estimates give only a rough idea of the computational 
demand of a method, we have looked also at the execution time, which is a simpler 
and a more informative assessment. 
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Method 

Complexity (single /) 

Complexity (vector of /) 

FFT2 

ff{N^\og{N)) 

ff{MN^log{N)) 

ENZ-PB 


^{N^M+N^KM + NK^^'^) 

EN-BB 

ff{N^K + NK^/^+K^) 

+ N^KM + NK^/^ + K^) 

ENZ-EBB 


ff{N^M + N'^KM+NK^/^ + 

GRBF 




Table 1: Estimates of the minimal computational complexity of the methods when evaluating U at 
a.N xN grid of nodes (r, 0) and for M values for the defocus parameter /, using K functions in the 
corresponding series expansions (for ENZ and GRBE). 


For the beginning, we run the ENZ-based and the GRBF algorithms evaluating 
17 at an 100 X 100 mesh of nodes for a single value of the parameter /, recording 
the execution time in dependence of the number of functions used in the corre¬ 
sponding series expansions. The experiments were performed on a standard PC 
running Matlab v. 8.1. Figure [T] shows the results, along with the corresponding 
regression curves. The slope of the regression line is approximately 0.0612 sec¬ 
onds/function for ENZ-PB, 0.0144 for ENZ-EBB and 0.0028 seconds/function for 
GRBF; although the dependence for ENZ-BB is clearly non-linear, the average 
slope of the regression curve in the studied range for this case is 0.2151. Thus, the 
GRBF approach is at least 5 times faster in comparison with its closest competitor, 
the ENZ-EBB method. For comparative purposes, the execution time to evaluate 
function U numerically making use of the two-dimensional fast Fourier transform 
was of approximately 0.25 seconds, matching the execution time for ENZ using 
about 12 Zemike terms, or for GRBF with approximately 90 Gaussian RBFs. 

In order to achieve a reasonable comparison, in the rest of our experiments we 
implemented a standard setting for each of the methods. Namely, the integral U 
was evaluated at a 100 x 100 regular grid of nodes by the semi-analytic methods. 
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Figure 1: Dependence of the computation time from the number of functions used in the descrip¬ 
tion of the complex pupil function, by means of the ENZ-PB (formulas (|34|l-l|35|l), the ENZ-BB 
(formulas (|41[l- (|4^), the ENZ-EBB (formulas 1(44^-1(47}) and the GRBE methods. 
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and at a 512 X 512 regular grid for the FFT2-approach (this is a realistic size needed 
to overcome aliasing and obtain accurate results). Moreover, we used K = Kgrbf = 
400 functions in the formula ( [24l ) (corresponding to the 20 x 20 regular grid of 
centers, as described above) and the cut-off parameter S = 60, while for the ENZ 
methods it was observed that a total of .S' = Kenz = 45 Zernike polynomials (up to 
the 8th order polynomials) was the optimal. This is actually about twice the number 
of terms recommended by ||8l HQ . Higher order polynomials, according to our 
experiments, did not contribute to higher accuracy, while causing ill conditioning 
of the computations. 

With these settings, we also compared the execution time of these methods as 
a function of the length of the vector of defocus parameters f. The results appear 
in Figure along with the regression lines for each scheme. The values of the 
slopes are approximately 0.139 for ENZ-PB, 0.064 for ENZ-BB, and for ENZ- 
EBB and 0.004 for GRBF, all in seconds per size of /. This means that ENZ-PB 
is about 35 times slower than GRBF when calculating U for many of values of f 
simultaneously, while ENZ-BB and EBB are about 16 times slower than GRBF, 
even though the number of Gaussian functions used was much higher than the 
number of Zemike polynomials. 

6.4. Accuracy 

The original ENZ theory, although representing a big step forward, had some 
limitations that had to be addressed. The obvious one is the use of only even terms 
in the Zernike expansion of the complex pupil function, which restricts it to the 
symmetric wavefront errors. Some other, less evident, problems lie in the core of 
the mathematical properties of the ENZ explicit formula for [/(r, 0;/). This is an 
infinite series of terms, each of them a finite linear combination of Bessel functions, 
and each new Zernike term added to the expansion of the pupil function ([Td]) in- 
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Length of defocus parameter vector 


Figure 2; Execution time according to the number of different defocus parameters used in the com¬ 
plex pupil function. A fixed number of functions (400 for GRBF and 45 for ENZ) was used for each 


value of /. 
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creases the complexity of the terms. The series is slowly convergent, especially for 
larger values of /, requiring a truncation with a large number of terms depending 
on / (it is recommended to use 3|/| + 5 terms, according to |[8ll26]l). 

Another issue in evaluating the ENZ expressions is the accuracy. The terms of 
the infinite series with even and odd orders form sign-changing sequences, which 
increases the risk of the cancellation errors. This phenomenon can be illustrated 
by the following experiments: in the case when the wavefront is given only by a 
positive Z| horizontal astigmatism, the evaluation of, say, the imaginary part of 
(7 at a point with radial coordinate r = 0.9 consists in adding a finite alternating 
sequence, with two dominant terms of approximately 0.423 each, with opposite 
signs, but whose absolute values differ in 3 x 10^^. This shows that these calcu¬ 
lations, if not well organized, can yield a loss of precision in about 5 significant 
digits. Last but not least, the ENZ formulas contain binomial numbers that must be 
evaluated with care in order to avoid overflow. 

Although the enhanced scheme of evaluation of the ENZ formulas ([44|)-(|47]) 
represents a significant improvement in speed, as the experiments show, still the 
danger of the cancellation errors is looming in the horizon due to the alternating 
sign terms in ( |44l ). 


In comparison, the error estimates from Section 4.3 show that the convergence 


in ( [2()| ) is very fast, and only a very reasonable number of terms is required for a 
precise evaluation. 

In order to assess accuracy, we carried out a number of numerical experiments 
where we calculated U. Eor the ideal wavefront (<I> = 0) and zero defocus (/ = 0), 
the results were discussed in llT/l . where it was shown that the two semi-analytic 
methods performed similarly, although GRBE calculations were done at a much 
lower cost. 

More informative is to use a synthetic wavefront described by a combination of 
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Zemike polynomial terms and exponential function, seeking a “fair” comparison 
between the ENZ (“Zernike-oriented”) and the GRBF (“exponentially-oriented”) 
approaches. Such class of functions allows also for an easy modeling of low and 
high-order oscillations. For experiments with a wavefront comprised of basically 
low-frequency oscillations see lITTll . Here we have used the function given by 

^{p,G) = 0.6Z|(p, e) - 0.4Z^(p, 0) - 0.3Z5 (p, 0) + 0.25Z^(p, 0) 

-hO.25Z^(p,0)-O.15Z^(p,0)-hO.4g(pcos0,psin0;-0.3,0,15) 

— 2 [g(pcos0,psin0;O.5,O.3,10) -|-g(pcos0,psin0;O.5, —0.3,10)], 

(49) 

where Z™ are the (orthonormal) Zemike polynomials and 

g{x,y;a,b,X) = exp{-X[{x-a)^ + {y-bf]}, 

as well as its “contaminated” version, where we have added to <I> a normally dis¬ 
tributed random noise of the form 0.5 * randn ( ) at a grid of 100 x 100 equally 
spaced points (see Figure [^. 

For these wavefronts we calculated the diffraction integral ( [TT] ) for different 
values of / by quadrature using the scientific software Mathematica with extended 
precision (with the options PrecisionGoal set to 8 and WorkingPrecision to 
16). These values of U (regarded as “exact”) were compared with the calculations 
performed by FFT2 and by two semi-analytic approaches discussed here, for which 
the diffraction integral was evaluated by all methods in a grid of 256 x 256 equally 
spaced points in the square [—2,2] x [—2,2]. 

Seeking a fair comparison of the computation of the integrals by both semi- 
analytic methods we fitted the pupil function using the first 200 Zemike polynomi¬ 
als for the ENZ, while for the GRBF the approximation was performed by the linear 
combination of 20 x 20 Gaussian functions, with the parameter A = 16. With this 
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Figure 3: 3D plot of the synthetic wavefront function defined in \‘X9) (upper picture) and of the same 
wavefront contaminated by a white noise 0.5 x N{Q, 1) (bottom). 
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choice of the number of functions in each basis an approximately similar goodness- 
of-fit to the wavefront data with noise is reached (the RMS error aproximately 0.09 
for both methods). 

As an illustration of comparative performance of the computation methods for 
different values of the defocus parameter / we plot in Figure the values of the 
normalized PSF, again for the simulated wavefront ( [4^ , along the horizontal line 
(0 = 0,71 and r G [0,1]), setting f = 0, f = In and / = —271, as in a numerical 
experiment described in f8]. In the case / = 0 the dotted line (computed by the 
ENZ method) is not observed because it matches exactly the other two curves, 
found using quadrature and the GRBF algorithm. We see that for larger values of 
defocus our procedure outperforms the alternative methods. 

In our last experiments we analyzed the performance of all methods for a non¬ 
circular geometry. Namely, for the same synthetic wavefront as before we set an 
elliptic pupil, taking in ( fT?] ) as A(p,0) the characteristic function of the set (in 
cartesian coordinates) 


x^ + - 


< 1 . 


(0.7)2 - 

The results are illustrated in Figure]^ The poorer performance of the ENZ method 
in these examples is due probably to a less accurate fit of the complex wavefront 
by Zernike polynomials over a non-circular domain. 


7, Conclusions 

A new procedure for computing 2D Eourier-type integrals, and in particular, the 
diffraction integrals with variable defocus has been put forward. We have discussed 
some error bounds and developed an efficient scheme for parallel evaluation of 
these integrals in a grid of points and for a vector of defocus parameters. 

It should be noted that when calculating the Eourier-type integrals considered 
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Section of PSF 



Figure 4: Values of the normalized PSFs and the corresponding errors (difference between the PSF 
computed by quadrature and by each other method, as explained in the text) for the wavefront \49\ 
+ noise, along the horizontal diameter of the unit disk, calculated for each method, for / = 0 (top 
row), f = 2tz (middle) and / = —2k (bottom row). 
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Figure 5: Values of the normalized PSFs for the synthetic wavefront with an elliptic pupil along the 
horizontal diameter of the unit disk, calculated for each method, for / = 0 (left) and f = 7l (right). 

in this paper by replacing the original integrand by its approximant, special care 
should be put in the quality of approximation. It is easily seen that proximity in 
an L? or even uniform norm is not enough to guarantee small errors (due to a high 
sensitivity of these integrals to oscillations), and Sobolev-type norms should be 
used. The higher accuracy we achieve in the approximation step, the better results 
will be obtained for the integral calculation. In our case, a linear least squares fit 
with Tikhonov regularization gave the best results. 

It is worth pointing out that the semi-analytic paradigm for computed Fourier- 
type integrals can be used with other bases, that might give other computational 
advantages, such as Sobolev orthogonal polynomials or polynomials introduced in 
ED, as well as for the computation with the more general focal factors, like in 

([T^. 

The proposed here GRBF-based approach has been compared with the two 
existing procedures, i.e., the 2-dimensional fast Fourier transform and the extended 
Nijboer-Zemike theory. The results of the comparison show that the new scheme is 
very competitive, providing higher accuracy and speed. Among other advantages 
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of our method we can mention the following: 


• a relative robustness of the computation with respect to the underlying ge¬ 
ometry of the pupil function. The GRBF have an almost local character 
(especially, for higher values of the shape parameter), and the fit that the lin¬ 
ear combination of such functions provides is less affected by the shape of 
the pupil; 

• the existence of the shape parameter in the GRBF provides more flexibility 
to small details, and allows for an easy implementation of a multi-resolution 
scheme, especially in the case of existence of subdomains with different 
complexity of the integrand. Since each function used for approximation 
of the pupil function enters the final expression linearly, one can use two 
or more layers of GRBF to fit the residual error consecutively using differ¬ 
ent sets of centers and different shape parameters in order to improve the 
accuracy of the results; 

• in the case of the calculation of diffraction integrals of the form ( [TT] ), the 
increase of the computational cost for a vector of values of the defocus pa¬ 
rameter is practically negligible, providing a substantial increase in the per¬ 
formance with respect to the other techniques. This is a reliable and efficient 
way of obtaining the through-focus characteristics of an optical system at 
higher resolutions in reasonable time. 

In conclusion, the proposed GRBF approach allows calculating through-focus 
point spread function at a very low computational cost for an arbitrarily selected 
set of the defocus parameters. This is particularly attractive in those applications in 
which evaluation of through-focus characteristics of an optical system is required. 
They include wavefront sensing, phase retrieval, lithography, microscopy, extreme 
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ultraviolet light optics, digital holography and physiological optics. We believe that 
the GRBF-based method has an even wider scope of application in mathematical 
imaging and vision. 
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